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ABSTRACT 

By analyzing numerically simulated solar oscillation data, we study the in- 
fluence of non-uniform distribution of acoustic wave amplitude, acoustic source 
strength, and perturbations of the sound speed on the shifts of acoustic travel 
times measured by the time-distance helioseismology method. It is found that 
for short distances, the contribution to the mean travel time shift caused by non- 
uniform distribution of acoustic sources in sunspots may be comparable to (but 
smaller than) the contribution from the sound speed perturbation in sunspots, 
and that it has the opposite sign to the sound-speed effect. This effect may cause 
some underestimation of the negative sound-speed perturbations in sunspots just 
below the surface, that was found in previous time-distance helioseismology in- 
ferences. This effect cannot be corrected by artificially increasing the amplitude 
of oscillations in sunspots. For large time-distance annuli, the non-uniform dis- 
tribution of wavefields does not have significant effects on the mean travel times, 
and thus the sound-speed inversion results. The measured travel time differ- 
ences, which are used to determine the mass flows beneath sunspots, can also be 
systematically shifted by this effect, but only in an insignificant magnitude. 

Subject headings: Sun: helioseismology - Sun: oscillations - sunspot 



1. Introduction 

Time-distance helioseismology is based on measuring and inverting acoustic wave travel 
times between separate points on the surface of the Sun. It is one of widely used approaches 
of local helioseismology for reconstructing solar subsurface structures and flows. Calculation 
of the temporal cross-covariance of two oscillation signals, observed in different points on 
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the so lar surface, is a key element of this method (IDuvall et al.lll993l ). iKosovichev fc Duvall 
( 119971 ) showed that the cross- covariance function for waves with the phase speed lying in 
a narrow interval can be approximately represented by a Gabor wavelet. The phase and 
group travel times of acoustic waves can be obtained by fitting the Gabor wavelet to the 
observed cross-covariance function, using a least-square method. The measured phase travel 
times are used for inferring the subphotospheric pertur bations of the sound (wave) speed 



and f low velocities in the quiet Sun and sunspot regions (IKosovichev et al.ll200d : IZhao et al. 



2OOII ). The mean travel time of acoustic waves traveling between two points in the opposite 
directions are used for determining the sound speed, and the travel time difference is used 
for determining the flows. 

However, an accurate inference of sunspot 's subsurface sound-speed structures and flow 
fields by use of this approach may be affected b y a series of phy sical and unphysical effects, 
such as strong wave damping in active re gions (IWoodardI 119971 ) . and the presence of strong 
magnetic field (IZhao &: KosovichevI |2006| ). Recently, by using observations of a quiet Sun 



region and artificially reducing solar acoustic oscillation amp litudes, i.e. masking solar wave- 
field to mimic the sunspot's behavior. iRajaguru et al.l (120061 ) found that this procedure could 



shift the measured acoustic travel times systematically by an amount of 5 — 40%, although 
such a shift was not expected. Furthermore, they suggested to correct the observed acoustic 
wavefields inside active regions by artificially increasing the wave amplitude. 

However, it is evident that the artificially masked wavefield of a solar quiet region 
can only mimic the acoustic power of the active region, but not the actual physical cause. 
Therefore, the systematic errors estimated by this approach may be inaccurate, and the 
correction procedure is unjustified. In this paper, we have carried out 3D numerical simula- 
tions of solar oscillations based on three different models to mimic the sunspot's wavefield, 
and investigated the systematic errors caused by the amplitude effects in the time-distance 
measurements. These r nodels include the artifi cial masking the numerically simulated wave- 
fields, as suggested by IRajaguru et al.l (120061 ) . reducing the stre ngth of oscillation sources 
to refl ect the physical effect of reduced excitation in sunspots (jParchevsky fc Kosovichev 
2007al ). and compare these with effects c aused by using a sound- speed perturbation deduced 
by the previous sound-speed inversions (IKosovichev et al.ll2000l ). The numerical simulation 
procedure and results are described in §2, and the results of the time-distance analysis are 
given in §3, followed by discussions in §4. 
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Numerical Simulations 



2.1. Governing Equations 



V ■ {pou) = 



Propagation of acoustic waves on the Sun is described by the system of the hnearized 
Euler equations 

d£ 
'dt 
d 

— (pow') + V p' = Oop' + f{x, y, z, t), 

where u' is the velocity perturbation, p' and p' are the density and pressure perturbations 
respectively, and f{x,y,z,t) is the function describing the acoustic sources. The pressure 
Pq, density po, and gravitational accelerations Qq with subscripts correspond to the back- 
ground model. To close the system ([1]) we used the adiabatic relation Sp/po = Ti6p/po 
between Lagrangian variations of pressure 6p and den sity 6p. The adiaba tic exponent Fi is 
calculated from the realistic OPAL equation of state (IRogers et al.lll996l ) for the hydrogen 
X and heavy elements Z abundan ces of the standard model. The standard solar model S 



( Christensen-Da. 



sgaard et al.lll996l ) with a smoothly joined model of the chromosphere of 



Vernazza et al.l (Il976l ) is used as the background model. 



The standard solar model is convectively unstable, especially in the superadiabatic sub- 
photospheric layers where convective motions are very intense and turbulent. Using this 
convectively unstable model as a background model leads to the instability of the solution 
of the linear system. The condition for stability against convection requires that the square 
of Brunt- Vaisala frequency Mq = go{l/Ti dlogpo/dr — dlogpo/dr) is positive. To make 
the background model convectively stable we replaced all negative values of //q by zeros 
and recalculated the profiles of pressure and density from the condition of hydrostatic equi- 
librium. This procedure guarantees convective stability of the background model. It has 
been shown that the profiles of pressure po, density po, sound speed cq, and acoustic cut-off 
frequency u;^ = Co/4if^(l — 2 dHp/dr) of the modified model are very close to the corre- 



sponding profiles of the standard solar model (jParchevsky fc Kosovichevll2007bl ). Quantity 
= —d\ogpo/dr represents the density hight scale. 

To prevent spurious reflections of acoustic waves from the boundaries we established non- 
reflecting boundary conditions based on the Perfectly Matched Layer (PML) method (IHu 
19961 ) at the top and bottom boundaries. The top boundary was set at the height of 500 km 
above the photosphere. This simulates a realistic situation when not all waves are reflected 
from the photosphere. Waves with frequencies higher than the acoustic cut-off frequency pass 
through the photosphere and are absorbed by the top boundary. This naturally introduces 
frequency dependence of the reflecting coefficient of the top boundary. The lateral boundary 
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conditions are periodic. The details are described by iParchevsky fc KosovichevI (j2007bl ). 
The waves are generated by spatially localized sources of the z-component of force 

.1 2 



f{x,y,z,t) 



with r and r given by 







'l-2T^)e-^ if r<R, 



if r > R, 



(2) 



^/{X- XsrcV + {y - VsrcY + - ZsrcY, T = ^ - VT, to < ^ < ^0 + — (3) 



UJ 



where is the unit vector in the vertical direction, Xsrc, Vsrc, and Zsrc are the coordinates of 
the center of the source, Rgrc is the source radius, u is the central frequency, to is the moment 
of the source ignition, A is the coefficient, which is measured in dyn cm~^ and describes the 
source strength. It has a physical meaning of the force density. 



2.2. Numerical Method 



To solve the system ([T]), a semi-discrete code developed by IParchevsky fc Kosovichev 



(l2007bl ) is used. The high-orde r dispersion relation preserving (DRP) finite difference (FD) 
scheme of iTam Sz WebbI fll993l ) is used for spatial discretization. The coefficients of this FD 
scheme are chosen from the requirement that the error in the Fourier transform of the spatial 
derivative is minimal. It can be shown that the 4th-order DRP FD scheme describes short 
waves more accurately than the classic 6th-order FD scheme. A 3rd-order, thre e-stage stro ng 
stability preserving Runge-Kutta scheme with the Courant number, C = 1, (jShul l2002l ) is 
used time advancing scheme. 

The efficiency of the high-order FD schemes can be reach ed only if they a r e corn bined 
with adequate numerical boundary conditions. We followed ICarpenter et al.l (119931 ) and 
used an implicit Pade approximation of the spatial derivatives near the top and bottom 
boundaries to derive a stable 3rd-order numerical boundary conditions consistent with the 
4th-order DRP numerical scheme for interior points of the computational domain. 

Waves with the wavelength less than 4Ax are not resolved by the FD scheme. They 
lead to point-to-point oscillations of the solution that can cause a numerical instability. Such 
waves have to be filtered out, and we use a 6th-order digital filter to eliminate unresolved 
short wave component from the solution at each time step. 
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2.3. Simulation of Artificial Data 

The simulations are carried out in a rectangular domain of size 90x90x31 Mm^ using 
a uniform 600x600x619 grid. The background model varies sharply in the region above 
the temperature minimum. Thus, to simulate the propagation of acoustic waves into the 
chromosphere we choose the vertical spatial step Az = 50 km in order to preserve the 
accuracy and numerical stability. The spatial intervals in the horizontal direction are Ax = 
Ay =150 km. To satisfy the Courant stability condition for the explicit scheme, the time 
step is set to be equal to 0.5 seconds. The sources of the z-component of force with random 
amplitudes and uniform frequency distribution in the range of 2 — 8 mHz are randomly 
distributed at the depth of 100 km below the photosphere and are independently excited at 
arbitrary moments of time. 

We describe three sets of simulations with different distribution of acoustic sources and 
different background models. The first reference model (Model I) represents simulations of 
the acoustic wave field for horizontally uniform distribution of the acoustic sources and the 
horizontally uniform background model. This model corresponds to the quiet Sun and will be 
used as a reference state for the following time-distance analysis. The acoustic travel times for 
models II and III are computed relatively to this reference model. The goal of this study is to 
estimate the contributions to the travel times arising from perturbations of the background 
model and non-uniform distribution of the acoustic sources separately. For this purpose, 
in model II the acoustic source strength is gradually decreased ( masked) in the central re- 



gion, s imulating the reduction of the acoustic sources in sunspots (jParchevsky fc Kosovichev 



2007al ). In this model, the horizontal axially symmetric distribution of the acoustic source 
strength is given by formula 

A(.,y)^{ K'^'^X) " (4) 

1 if r/j > Rs, 



where = -y/(x"^^a?^p~+7y"^-ys)^ is the horizontal distance from the sunspot axis, Xs, Us, 
and Rs are the x-, y- coordinates of the sunspot center and the sunspot radius respectively. 
The background model remains unperturbed and horizontally uniform. So, as far as the 
background model remains unchanged all deviations of simulated wave field properties from 
model I can be explained as a result on non-uniform distribution of the acoustic sources. 

Strictly speaking, travel times for the cases of uniform distribution of the acoustic sources 
(Model I) and masked source strength (Model II) are calculated at different conditions. The 
amplitude of the wave field is uniform in the first case and non-uniform in the second one. To 
take this into account we mask the wave field of the Model I by masking function computed by 
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averaging signals azimuthally around the su nspot center. Th i s mirn ics the reduced amphtude 
of active regions, just as what was done by lRajaguru et all (120061 ) . The resultant model is 
called Model la for the convenience of reference in the following descriptions. Although 
the amplitude distributions now are the same for both wave fields, the wave fields itself are 
different, because masking of the source strength is not reduced to simple masking of the 
resulting wave field. 

Model III combines the source masking of Model II with the sound speed perturbation 
in sunspots. The 3D sound-speed profile c in Model III is approximated by the formula 



cix,y,z) = co{z) 



6c{z) 
coiz) 



Aix,y)) 



(5) 



where 6c/ cq is the vertical profile of the sound speed perturbation at the sunspot axis. This 
profile is shown in Figure [Hand was calculated from the inversion of helioseis mic data for the 
sunspot observed by SOHO/MDI on 20 June 1998 flKosovichev et al.ll2000l ). The change of 
the sound-speed perturbation sign from negative to positive at approximately 4 Mm below 
the photosphere is a characteristic feature of this profile. The depth of inversion divides the 
domain in vertical direction into two regions with the sound speed greater and smaller than 
in the standard reference model. Hence, we expect different behavior of waves propagating 
though this artificial sunspot if their turning points lie in different regions. 

The amplitude map of the resulting wave field for Model III is shown in the left panel 
of Figure [21 The solid line in the right panel represents the azimuthally averaged amplitude 
profile. The dashed line shows the angularly averaged amplitude profile for Model II. The 
inhomogeneity of the sound speed causes increasing the ratio of oscillation amplitudes outside 
and inside of the artificial sunspot by about 40%. 

The acoustic power spectrum {k-u diagram) of the simulated wave field is shown in 
Figure [31 We see a good agreement with the observed power spectrum in terms of the 
shape and locations of power ridges, yet the simulated wave field has more power in high- 
frequency region. It is not clear whether such power excess really exists and instrumentally 
filtered out during observations or it is an artifact of the numerical modeling of the wave 
field. The realistic non-linear simulati ons of solar convection also show a power excess in 
the k-u diagram at higher frequencies (iGeorgobiani et al.l 120071 ). For the present study, the 
high-frequency power excess is unimportant, because for the time-distance helioseismology 
analysis, we only select the frequency band of 3 - 4 mHz for analysis, and in this range, the 
simulated acoustic power is similar to the observations. 
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3. Results 
3.1. Time-Distance Analysis 

To perform the time-d istance helioseismo losy analysis of the simulated data we followed 



the procedure described in IZhao et al.l (120011 ). For all our models, we select the annuli with 
radii of 6.2 — 11.2 Mm, 8.7— 14.5 Mm, and 14.5 — 19.4 Mm to obtain the mean travel times, 
Tmean, and the travel time differences, r^^g, (index i marks model number, including the 
reference Model 1), which are respectively averages and differences of outgoing and ingoing 
travel times in the time-distance center-annulus measurement scheme. A phase-speed filter is 
applied in each case to select only waves in a narrow phase speed interval. To study the effects 
caused by wavefield non-uniformity, we calculate the differences of 6Tmla.n = TmLn — TmLn 
and 5r^ig = r^-g — r^|^ for the analysis. 

In Figure m we show the maps of the mean travel time perturbations, (5rmean, the travel 
time differences, Sr^^g, for all three models. Although the background model of Models la 
and II are the same with the reference Model I, we can see systematic shifts of mean travel 
times inside the masked regions. 

For better understanding the results, it is useful to compare the profiles of the travel 
time deviations, azimuthally averaged around the sunspot center for both ^TmLn and St^^q, 
as shown in the middle row of Figure [51 Obviously, the mean travel time shifts, (5rmean, 
are significantly larger in Model II than in Model la, although both have exactly the same 
background model and exactly the same oscillation amplitude reduction in wavefields. Ex- 
pectedly. Model III shows mostly positive travel time shifts in contrast with the other two 
models, and this is obviously due to the negative sound-speed perturbation to the back- 
ground model in a shallow subsurface region. One would expect this positive time shift 
would increase significantly if there is no such effect that causes the time deficit in Model II, 
however, it is not immediately clear whether this is the case, or if it is, how much it would 
increase. For the travel-time difference, STdm, the shifts for all three models are quite small, 
within an order of 2 sec, substantia lly smaller than the measured travel time shifts from a 
real sunspot data ( Zhao et"aDl200ll ). 



The azimuthally averaged (5rmean and ^t^-^ for the other two annuli measurements are 
also presented in Figure [51 For the shorter travel distances, both Models la and II show 
stronger travel time deficits in the mean travel time measurements compared the the inter- 
mediate travel distance case, up to approximately 15 sec for Model II. However, Model III 
still displays mostly a positive sign, although it displays some dips in the central area where 
one would expect a stronger positive shift because of the larger negative sound-speed pertur- 
bation there. Again, Srj^^ does not show any significant time shifts for all models. For the 
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larger annulus radius, Models la and II do not display significant time shifts, but Model III 
displays a significant negative time shift because for this set of measurement, waves reach 
the depth of a large positive sound-speed perturbation. For this annulus measurement. Mod- 
els II and III display an order of magnitude of 5 sec travel time shifts in ^Tdiff, larger than 
those of shorter annuli measurements, but still significantly smaller than the shifts in the 
real sunspot measurements. 



3.2. Power Correction 

Based on their artificial tests with the quiet Sun data, in or der to remove the mea sured 



travel time shifts caused by the oscillation amplitude reductions, iRajaguru et al.1 (120061 ) have 
suggested to make corrections for these areas by enlarging the observed oscillation amplitude 
in active regions. This procedure is just reverse to the artificial masking. It obviously works 
if the power reduction is caused by artificial masking, like in Model la. However, it is useful 
to examine whether this works for the oscillation power reductions that are not caused by 
surface masking, but physical mechanisms, such as the reduction in the excitation power 
(like in Models II and III). 

For each model (la, II, and III) we calculate the average amplitude profile and normal- 
ized the wavefield by using this profile (procedure of unmasking the wavefield), making the 
oscillation power nearly uniform over the whole box. The same time-distance analysis is 
performed like in §3.1, and azimuthally averaged curves are displayed in Figure |6l It can be 
clearly seen that, as expected, this power correction removes all travel time shifts in both 
^Tmean cLud ^idiff for all anuuli measurements for Model la. For ^Tmean, for the two shorter 
annuli measurements, the correction slightly lifts both Models II and III without changing 
signs of the profiles, and for the longest annulus, the correction does not change much the 
measurements. For 6Tdis, the correction changes the profiles of Models II and III for all 
annuli, but still, the travel time shifts are within 5 sec or so. 



4. Discussion and Conclusion 

The explanation of the fact that the acoustic travel times depend on the non-uniformness 
of the wave field amplitude or non- uniform distribution of the source strength is related to the 
definition of travel times in helioseismology, which have to deal with the stochastic randomly 
excited oscillations, rather than isolated point sources. The travel time of a wave packet 
traveling between two points on the surface is defined not as a local physical quantity, which 
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can be explicitly computed from the background model, but rather in an "observational" way 
as a parameter which is obtained from fitting of cross-correlation of the oscillation signals 
by the Gabor wavelet. Thus, it is very important to investigate the effect of non- uniform 
distribution of acoustic sources, damping and other causes of the non-uniform wave field 
distribution on the Sun. We have presented the results for some of these effects by using 
numerical 3D simulations of acoustic wave propagation in various solar models. 

We have found that the source masking for horizontally uniform background model 
(dashed curves) may cause a systematic negative shift of about 8-13 seconds in the mean 
travel times for short distances (annuli with radii smaller than 14.5 Mm). Such travel time 
shift may cause underestimation of the sound speed perturbation in the shallow (1-2 Mm 
deep) subsurface layers. For larger distances, the contribution to the mean travel time shift 
becomes negligible. On the contrary, the shift of the travel time differences (due to the 
non-uniform distribution of the acoustic sources) is negligible for short distances and has a 
value about —5 sec for the largest distance used in our experiments. This is much smaller 
than perturbations of the travel-time differences observed in real sunspots. 



The results of our experiments are different from a similar work by iHanasoge et al. 



( 120071 ). where authors report significant disbalance between ingoing and outgoing travel times 
(about —5 s for distance of 6.2 Mm and about —15 s for distance of 24.35 Mm), and suggested 
that at the large distances the false travel time difference signal caused by non-uniform 
distribution of sources may be misinterpreted as a result of subsurface flows. However, it is 
quite clear that, as shown in the first two annuli measurements of Figure the oscillation 
power deficit due to the source masking may have greatly reduced the travel time shifts 
measured in Model III, which means that if doing inversions, the inverted sound-speed profile 
would be gr eatly underestimated. T his suggests that the sound-speed profile under sunspots 



obtained by iKosovichev et al.l (120001 ) got the correct sign but might be underestimated. For 
the flow fields, this masking effect might cause some systematic velocity errors, but only of 
a very small magnitude. 

In addition, our experiments show that the amplitude reduction is caused by the weaker 
oscillation sources in sunspots cannot be corrected by a simple normalization procedure. This 
imposes us a difficult task on how to retrieve accurately the sound-speed profiles beneath 
sunspots, and improve the time-distance helioseismology inferences. 



REFERENCES 

Carpenter, M.H., Gottlieb, D., and Abarbanel, S. 1993, J. Comput. Phys. 108, 272. 



- 10 - 

Christensen-Dalsgaard, J., Dappen, W., Ajukov, S. V., Anderson, E. R., et al. 1996, Science, 
272, 1286. 

Duvall, T. L., Jr., Jefferies, S. M., Harvey, J. W., & Pomerantz, M. A. 1993, Nature, 362, 
430. 

Georgobiani, D., Zhao, J., Kosovichev, A. G., Benson, D., Stein, R. F., & Nordlund, A. 
2007, ApJ, 657, 1157. 

Hanasoge, S. M., Couvidat, S., & Rajaguru, S. P. 2007, ApJ, submitted 

Hu, F. Q. 1996, J. Comput. Phys., 129, 201. 

Kosovichev, A. G. & Duvall, T. L., Jr. 1997 in Proc. SCORe'96: Solar Convection and 
Oscillations and their Relationship, eds. F. P. Pijpers, J. Christensen-Dalsgaard, and 
C. S. Rosenthal, Kluwer Academic Publishers, Dordrecht, Holland, 241. 

Kosovichev, A. G., Duvall, T. L. Jr., & Sherrer, P. H. 2000, Solar Phys. 192, 159. 

Parchevsky, K. V., & Kosovichev, A. G. 2007a, ApJ, 666, L53 

Parchevsky, K. V. & Kosovichev, A. G. 2007b, ApJ, 666, 547. 

Rajaguru, S. P., Birch, A. C, Duvall, T. L. Jr., Thompson, M. J., & Zhao, J. 2006, ApJ, 
646, 543. 

Rogers, F. J., Swenson, F. J., & Iglesias, C. A. 1996, ApJ, 456, 902. 

Shu, C.-W. 2002, in Collected Lectures on the Preservation of Stability under Discretization, 
eds. D. Estep and S. Tavener, SIAM, 51. 

Tam, C. and Webb, J. 1993, J. Comput. Phys. 107, 262 

Vernazza, J. E., Avrett, E. H., and Loeser, R. 1976, ApJS, 30, 1. 

Woodard, M. F. 1997, ApJ, 485, 890 

Zhao, J., & Kosovichev, A. G. 2006, ApJ, 643, 1317 

Zhao, J., Kosovichev, A. G., & Duvall, T. L.,Jr. 2001, ApJ, 557, 384. 



This preprint was prepared with the AAS I^TJt;X macros v5.2. 



- 11 - 




Fig. 1. — Perturbation of the sound speed profile at tlie sunspot axis, calculated by inverting 
data obtained from SOHO/MDI observations of the sunspot on 20 June 1998. 
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K, Tilm X, Mm 

Fig. 2. — Amplitude map {left) of the vertical velocity component for simulations with 
masked sources and sound speed perturbation of the background model. The solid line 
represents azimuthally averaged profile of the map. The dashed line shows azimuthally 
averaged profile of the vertical velocity component for simulations with masked sources only 
(the background model is horizontally uniform). 
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Fig. 3. — Power spectrum diagrams from MDI high resolution observations (left) and Model I 
of our simulations (right). Both power spectra are computed using the same time duration 
and after the simulated data are binned down to the same spatial resolution as the observed 
data. 
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Fig. 4. — Maps for mean travel times {upper row) and travel time differences {bottom row) 
relative to the Model I. For both rows, from the left column to the right are for Model la, 
Model II, and Model III, respectively. Annulus radii used for time-distance measurements 
are from 8.7 to 14.5 Mm. 
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Fig. 5. — Mean travel times {left panels) and travel time differences {right panels) azimuthally 
averaged from maps like shown in Figure H] for different annulus radii of (from above to 
bottom) 6.2 - 11.2 Mm, 8.7 - 14.5 Mm, and 14.5 - 19.4 Mm. 
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Same as Figure [5], but after the power corrections. 



